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O ■ Abstract 



We present algorithms for computing ranks and order statistics in the Farey sequence, taking time 
0(n 2/3 ). This improves on the recent algorithms of Pawlewicz [Paw07| . running in time 0(n 3/4 ). We 
also initiate the study of a more general algorithmic problem: counting primitive lattice points in planar 
shapes. 



Since the publication of this technical report, this work has been extended and merged with 
the paper of Pawlewicz. The merged version is available at: 

■ http://web.mit.edu/~mip/www/papers/farey2/paper.pdf 

1 An Improved Algorithm for the Farey Sequence 

£Nj ■ The Farey sequence of order n, denoted !F n , is the ordered list of irreducible fractions £ with a < b <n. This 

sequence is a well-studied mathematical object, with fascinating properties. Sec GKP94] for a discussion at 
length. The sequence has 0(n 2 ) terms, and there are many algorithms for generating it entirely in 0(n 2 ) 
time. Perhaps the best known ones are based on the Stern-Brocot tree, and the properties of the mediant. 
. One can ask, however, for more local access to the sequence. Two natural questions arise: 



given a number x £ [0, 1] find Rank(x, n) = \T n n [0, x] 



OO 

o 

• given an index k < |.F n | find Statistic^, n) = the fc-th value in T n (in sorted order). 



The Statistic problem can be solved with O(lgn) calls to the Rank problem |Paw07j . so below only 
bounds for the Rank version are discussed. 

To the best of my knowledge, the question was first formulated in 2003, when I proposed it as a contest 
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problem at the 11th Balkan Olympiad in Informatics. The official solution consisted of an 0(n lg n) algorithm, 
and several contestant also found this algorithm. We |PP04j later described a solution with a slightly better 
running time of 0(n). 

Quite recently, Pawlewicz [Paw07] broke the linear time barrier, and provided an algorithm with running 
time 0(n 3 / 4 ). In the present, I describe an improved 0(n 2 / 3 lg 1 ^ 3 n) algorithm. 

1.1 Review: The Algorithm of Pawlewicz 

Let S n (x) = \{f \ b < n A | < x A gcd(a, b) — l}|. This provides the quantity we want to compute. It 
can be seen that: 

n 

S n (x)=Y\bx\-Y J S V ^ x ) (!) 

h=l d>2 

It is shown in |Paw07| that A n (x) = 1 bx \ can be computed in O(lgn) time. 
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Thus, the only challenge is to estimate the recursive component of the sum. The recursion will only need 
S ln/dl (x) for all d, since [[% JAfeJ = LdfkJ ■ 

The crux of the algorithm is the following observation: given all relevant Si(x), for i < k, then Sk(x) can 
be computed in 0{\fk) time. Indeed, J2d>2 S[k/d] ( x ) only contains at most 1\fk distinct terms: \fk terms 
corresponding to d < yk, and at most \fk terms for d > \fk because then we have k/d < \fk. The latter 
terms have multiplicities, but the multiplicity of each term can be computed easily in 0(1) time. 

Applying this observation to compute all needed S terms recursively, the running time is: 

n i — \fn J — </n </n ^/n 

d=l » d=l " 3=1 d=l * 3=1 

1.2 Our Improved Algorithm 

The key to our improved algorithm is to show that Si, . . . , Sk can be computed in time 0{k lg k). Then, the 
remaining terms can be computed by the old algorithm in time \J n l & = V™' ' 0(y/n/k) = 0(n/Vk). 

The total running time is then 0(fclgfc + n/yk), so it is optimized by picking k = (n/lgn) 2 / 3 . Thus, the 
running time is 0(n 2 / 3 lg 1 ^ 3 n). 

To compute Si, . . . , Sk efficiently, we make the observation that the composition (with respect to recursive 
terms) of Si and Si-i are not too different. Specifically, J2d>2 S[i/d\ ( x ) differs from J2d>2 S[(i-i)/d\ (x) only 
for the values of d that i is a multiple of. 

To maintain understanding of divisibility by all d's, and compute Si(x) values in order, we use an 
algorithm similar in flavor to Eratosthene's sieve |Hor72] . We first create an array D[l . . k], where D[i] holds 
a list of all divisors of i. To create the array, simple consider all d, and add d to D[md], for all m. This takes 
timeE^^Offclgfc). 

Now, iterate i from 1 to fc, maintaining J2d>2 ^N/ d J a ^ a ^ ^i mes - The sum is updated by considering 
all d € D[i], subtracting S\a_iyd\ (x) and adding Su/m (x). The complexity is linear in the size of D[l . . k], 
and is thus O(fclgfc). To compute Si(x) from this running sum, we only need Ai(x), which takes O(lgi), 
giving an additive O(fclgfc). 

2 Counting Primitive Lattice Points 

A primitive lattice point is a point (x,y) in the plane, with x,y € Z and gcd(a;,y) = 1. We observe that 
computing Rank (a;, n) in the Farey sequence is equivalent to couting primitive lattice points inside the right 
triangle defined by (0, 0), (n, 0) and (n, xn). Let us now generalize the algorithm to counting primitive lattice 
points in more general shapes. 

Counting primitive lattice points in planar shape is a relatively new topic in mathematics, but one that 
is gathering si gnifican t momentum |Mor851 INow881 IHen941 INow95[ IHN961 IMul96L INow971 IKN011 IZC991 
BCZ00 ( Wu02, Zha03, Now05 . In the mathematical sense, "counting" refers to estimating the number of 
primitive points with a small error, as the size of the shape goes to infinity. 

In this paper, we initiate the study of the algorithmic problem of counting (exactly) the number of 
primitive lattice points inside a given shape. More precisely, we study this problem for polygons containing 
the origin. The condition that the shape should contain the origin also appears in the mathematical works 
referenced above, and is natural given that we are counting points visible from the origin. 

Theorem 1. Let P be a polygon containing the origin, defined by k vertices at b-bit rational coordinates. 
If D is the diameter of the polygon, one can count the number of primitive lattice points inside P in time 

//• ■ - k-ir . 

A very pertinent question is how efficient this running time actually is. Remember that |PP04j shows 
that the Farey rank problem can be used to factor integers. Since counting primitive lattice points is a 
generalization, we conclude that a polynomial-time algorithm, i.e. poly (A; -6), is likely impossible. Thus, 
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the algorithm needs to depend on some parameter describing the polygon, which can be exponential in b. 
One such parameter is the diameter D. Clearly, however, this is not the only choice, and it is conceivable 
that other measures lead to better results. For right triangles such as those in the Farey rank problem, the 
diameter is n, yet we know a better algorithm with time essentially n 2 / 3 . 

A trivial alternative to Theorem[T]is the algorithm which iterates over all lattice points inside the polygon, 
and runs Euclid's algorithm on each point. It is possible [Bar94. MB02. Hir05] to list all lattice points with 
a polynomial poly(fe • b) cost per point. If the polygon has integral coordinates, Pick's formula shows that 
the number of lattice points inside is asymptotically equal to the area. Thus, the exhaustive algorithm has 
complexity proportional to the area, times polynomial factors. 

Unfortunately, the area and the diameter are not related in the worst-case (e.g., for very skinny shapes). 
However, in the more "typical" case when the polygon is fat, the area is A = Q(D 2 ). Thus our running time 
of D 6 / 7 can be rewritten as A 3 / 7 , which gives a significant saving over the exhaustive algorithm. It is an 
interesting open problem to construct an algorithm which beats exhaustive search for any polygon. 

A standard idea. Let P be a polygon, defined by rational points (x±,y±), . . . , (xk, yk)- Then, let Pu be 
the polygon defined by (^f, ^-), . . . , Define A(P) to be the number of lattice points inside polygon 

P, and S(P) the number of primitive lattice points inside P. 
We first observe the following recursive formula: 

S(P)=A(P)-J2S(P/ d ) (2) 

d>2 

Indeed, every point in A(P), but not in S(P) is a nonprimitive lattice point (x,y). If gcd(x,y) = d > 1, 
then (t, K) is a primitive lattice point. Furthermore, such a point is inside Pu, so we can remove all points 
with a greater common divisor equal to d by subtracting Pf d . 

We can bound the recursion depth in by appealing to the diameter D. We note that the diameter of 
P/(D+i) is l ess than 1, so it does not contain any lattice point outside the origin. Thus, S(P/^) = for all 
D < d < oo, and S(P/ 00 ) = 1 (the origin). Then, it suffices to consider only P, P/ 2 , ■ ■ ■ , P/d i n the algorithm. 

To compute the "constants" A{Pn) in the recursion, one needs to compute the number of lattice points 
inside polygons with rational coordinates. This is a well-studied problem, and [Bar94| MB021 IHir05j given 
polynomial-time algorithms. Observe that for i < D < 2 b , coordinates of Pn have at most 2b bits of 
precision. Thus, computing A(Pu) takes time 0(k ■ poly(fr)); the dependence on k is linear by triangulating 
the polygon. 

Note that formula @ is very similar to the recursive formula for the Farey rank problem JT]). Indeed, 
is simply a transcription of (HJ), where the polygons are the relevant right triangles. 
It is thus tempting to conjecture that we can use the same dynamic program to evaluate ([2]). Unfortu- 
nately, this is not the case, for a somewhat subtle reason. Due to the geometry of the rank problem, SL/<j 
was the same as S\ n u\ . By taking floors, we concluded that among S\ n u\ 's with d 6 {^/ri, ■ ■ ■ , n}, there are 
only y/n distinct quantities. Unfortunately, in general we cannot round the vertices of P14 to lattice points, 
and thus we cannot conclude that for d > \f~D there are only \f~D distinct cases. 

A more careful analysis. Since we are not interested in factors of k ■ b°^ in the running time, let us 
define 0*(f) = f ■ k ■ . To use ideas similar to the previous dynamic program, we begin by obtaining a 
weaker bound for computing the small terms of the recurrence: 

Lemma 2. We can (implicitly) compute S(P/ T ), S(P// T +i\), S(Pm) in time O* (t-) ■ 

Proof. Note that P/d Q P/(d-i) Q ' ' ' Q P/t- Also, the diameter of Pi r is D/t, implying that P/ T , and any 
smaller polygon only contain 0((r=r) 2 ) lattice points. Since S(P/i) can only grow between 1 and 0((-^) 2 ) 
when i goes from D to r, there are only so many distinct values that can appear. 

To actually compute these values, we perform an exhaustive enumeration of all lattice points in P/ T . 
Points which are not primitive are discarded. For every primitive point (x,y), we compute a value ip(x,y) 
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which is the minimum i such that Pu does not contain it. This is done by a binary search for i. Every 
comparison is a point-in-polygon test, which takes 0{k) time. 

We now sort the (j)(x, y) values from largest to smallest. The sorted list gives us an implicit representation 
for S(P/i) for every i > r. Indeed, it suffices to binary search for the first occurrence of i; the elements to 
the left correspond to primitive points which are inside P^. □ 

Let us now consider the problem of computing a term using our recursion: S(Pu) = A{P/ i )—^2 d>2 S{P/ i( i). 
We assume terms S{P/j) for j > i have already been computed; in particular, for j > r we have the implicit 
representation from Lemma [2] 

In principle, the expression for S(P/i) has \Jj\ terms. However, as in Lemma [U we can observe that for 

fixed 6, there are only O(-^r) distinct terms among all S(P/ig)'s with id > S. These terms can actually be 
summed up in O* (-jt) time. Indeed, we begin with do = |~|] , and binary search for the minimum d\ such 
that S(P/id 1 ) < S(P/id a )- We add to the sum S(P/id ) ■ (di — do), then binary search for the minimum di 
leading to a different value, and so on. 

After dealing like this with all terms d > 4, we can simply add the remaining terms. The running time 

for computing S(P/i) is therefore O* (-p- + |). We can optimize by setting S = D 2 /^ 1 / 3 , yielding a running 
time of O* ((f) 2 / 3 ). 

We now want to optimize the parameter r in Lemma [2] A choice of r implies that we spend 0*(^r) 
time by Lemma [5J and then run the dynamic program for terms S(Pu) with i < r. This second part will 
take time: 

/ n\ 2/3 

Then, we can optimize by setting = D 2 I 3 t x I 3 , i.e. r = D 4 / 7 . The final running time is therefore 0*{D e / 7 ). 
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